function [S, dS] = predict_IRT1(theta, t)
% [S, dS] = predict_IRT1(theta, t)
% Predicts the signal intensities of magnitude MR images in a
% inversion recovery T1 measurement.
%
% Model:
%   S(t) = | A - B exp(- t * (R1/1000)) |
% 
%   R1 = theta(1)
%   A = theta(2) 
%   B = theta(3)
%
% Created by Dirk Poot Erasmus MC, 
% 26-8-2011

R1 = theta(1,:)/1000;
A = theta(2,:);
if size(theta,1)>2
    B = theta(3,:);
else
    B = 2*A;
end;

% S = A -B * exp(-t*R1)
E1 = exp(-t*R1);
S = bsxfun(@plus, A , bsxfun(@times, - B , E1));

if nargout>1
    sgn = sign(S);
    % also compute derivative:
    dSdR1 = (t*B) .*E1 .* sgn;
    dSdA  = sgn;
    dSdB  = -E1.*sgn;
    if size(theta,1)>2
        dS = cat(3,dSdR1/1000, dSdA, dSdB);
    else
        dS = cat(3,dSdR1/1000, dSdA+2*dSdB);
    end;
end;
S = abs(S);